PT_bias
Fernando Rawicz
12/05/2021
options(scipen = 9999)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
set.seed(415)
library(rstan)
library(shinystan)
library(stringr)
library(dplyr)
library(readxl)
library(gtools)
library(rjags)
library(runjags)
library(brms)
library(knitr)
library(DT)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())R Markdown
candidates <- c("PS", "CHEGA", "PCP", "PSD", "BE", "IL", "RIR")
pt_pitagorica <- read_excel("data/pollsportugal.xltx", sheet = "Pitagorica")%>%
filter(Party %in% candidates)%>%
mutate(Bias = Result - Poll ) %>%
mutate(Company = "Pitagorica")
pt_UCP <- read_excel("data/pollsportugal.xltx", sheet = "UCP")%>%
filter(Party %in% candidates)%>%
mutate(Bias = Result - Poll ) %>%
mutate(Company = "UCP")
pt_eurosondagem <- read_excel("data/pollsportugal.xltx", sheet = "Eurosondagem")%>%
filter(Party %in% candidates)%>%
mutate(Bias = Result - Poll ) %>%
mutate(Company = "Eurosondagem")
pt_intercampus <- read_excel("data/pollsportugal.xltx", sheet = "Intercampus")%>%
filter(Party %in% candidates)%>%
mutate(Bias = Result - Poll ) %>%
mutate(Company = "Intercampus")
pt_aximage <- read_excel("data/pollsportugal.xltx", sheet = "Aximage")%>%
filter(Party %in% candidates)%>%
mutate(Bias = Result - Poll ) %>%
mutate(Company = "Aximage")
pt_bias_final <- rbind(pt_pitagorica, pt_UCP, pt_eurosondagem, pt_intercampus, pt_aximage) %>%
mutate(Companyiid = match(Company , sort(unique(Company)))) %>%
mutate(Eleciid = match(Election , sort(unique(Election)))) %>%
mutate(Partyiid = match(Party , sort(unique(Party))))
DT::datatable(pt_bias_final)boxplot(Bias ~ Company, data=pt_bias_final, main="Boxplot Bias all elections portuguese companies")
boxplot(Bias ~ Company, data=pt_bias_final %>% filter(Election == "Pres"), main="Boxplot Bias President portuguese companies")
Model with just house and party
mod_string = "
model {
## sampling
for (i in 1:N){
y[i] ~ dnorm(mu_party[PartyIndex[i]] + mu_house[HouseIndex[i]] , invsigma2)
}
## priors
for (j in 1:J){
mu_party[j] ~ dnorm(mu, invtau2)
}
for (h in 1:H){
mu_house[h] ~ dnorm(mu, invtau2)
}
invsigma2 ~ dgamma(1.5, 6)
sigma <- sqrt(pow(invsigma2, -1))
## hyperpriors
mu ~ dnorm(0, 4)
invtau2 ~ dgamma(3, 1)
tau <- sqrt(pow(invtau2, -1))
}"
params = c("mu_party", "mu_house", "mu", "invtau2", "invsigma2")
the_data <- list("y" = pt_bias_final$Bias,
"PartyIndex" = pt_bias_final$Partyiid,
"HouseIndex" = pt_bias_final$Companyiid,
"N" = length(pt_bias_final$Bias),
"J" = length(unique(pt_bias_final$Partyiid)),
"H" = length(unique(pt_bias_final$Companyiid))
)params = c("mu_party", "mu_house", "mu", "invtau2", "invsigma2")
mod = jags.model(textConnection(mod_string), data=the_data, n.chains=3)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 69
## Unobserved stochastic nodes: 14
## Total graph size: 260
##
## Initializing modelupdate(mod, 1e3)
mod_sim = coda.samples(model=mod,
variable.names=params,
n.iter=5e3)
mod_csim = as.mcmc(do.call(rbind, mod_sim))
plot(mod_sim)
traceplot(mod_sim, col = 1:3)
gelman.diag(mod_sim)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## invsigma2 1 1
## invtau2 1 1
## mu 1 1
## mu_house[1] 1 1
## mu_house[2] 1 1
## mu_house[3] 1 1
## mu_house[4] 1 1
## mu_house[5] 1 1
## mu_party[1] 1 1
## mu_party[2] 1 1
## mu_party[3] 1 1
## mu_party[4] 1 1
## mu_party[5] 1 1
## mu_party[6] 1 1
##
## Multivariate psrf
##
## 1dic = dic.samples(mod, n.iter=1e3)
(pm_params = colMeans(mod_csim))## invsigma2 invtau2 mu mu_house[1] mu_house[2] mu_house[3]
## 0.18529063 2.83215636 -0.28701191 -0.22371688 -0.34520823 -0.29310483
## mu_house[4] mu_house[5] mu_party[1] mu_party[2] mu_party[3] mu_party[4]
## 0.09527954 -0.91251443 0.16557152 -0.18953359 -0.12248021 0.05642923
## mu_party[5] mu_party[6]
## -1.07897684 -0.79943599summary(mod_sim)##
## Iterations = 1001:6000
## Thinning interval = 1
## Number of chains = 3
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## invsigma2 0.18529 0.0325 0.0002654 0.000294
## invtau2 2.83216 1.4721 0.0120197 0.021036
## mu -0.28701 0.2411 0.0019683 0.003045
## mu_house[1] -0.22372 0.5045 0.0041192 0.004898
## mu_house[2] -0.34521 0.4972 0.0040597 0.004686
## mu_house[3] -0.29310 0.5229 0.0042698 0.004971
## mu_house[4] 0.09528 0.5499 0.0044898 0.005681
## mu_house[5] -0.91251 0.4113 0.0033586 0.004418
## mu_party[1] 0.16557 0.4729 0.0038609 0.005071
## mu_party[2] -0.18953 0.6314 0.0051553 0.005968
## mu_party[3] -0.12248 0.6261 0.0051121 0.005772
## mu_party[4] 0.05643 0.6260 0.0051116 0.006400
## mu_party[5] -1.07898 0.4620 0.0037724 0.004992
## mu_party[6] -0.79944 0.4364 0.0035635 0.004503
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## invsigma2 0.1281 0.1623 0.18338 0.20583 0.25516
## invtau2 0.8791 1.7814 2.53183 3.55290 6.48869
## mu -0.7461 -0.4490 -0.29317 -0.12875 0.20688
## mu_house[1] -1.1972 -0.5587 -0.23010 0.10560 0.79556
## mu_house[2] -1.3202 -0.6709 -0.35182 -0.02563 0.66141
## mu_house[3] -1.3200 -0.6375 -0.29779 0.04362 0.76527
## mu_house[4] -0.9287 -0.2806 0.07539 0.44576 1.22297
## mu_house[5] -1.7501 -1.1803 -0.89741 -0.63392 -0.12991
## mu_party[1] -0.7158 -0.1585 0.15043 0.47293 1.12995
## mu_party[2] -1.4283 -0.6053 -0.20442 0.20857 1.09561
## mu_party[3] -1.3096 -0.5367 -0.14784 0.27320 1.18168
## mu_party[4] -1.0920 -0.3611 0.02871 0.43778 1.37623
## mu_party[5] -2.0235 -1.3776 -1.06773 -0.76618 -0.20058
## mu_party[6] -1.6898 -1.0789 -0.78776 -0.50435 0.02468aaa <- summary(mod_sim)
# xtable(as.data.frame(aaa$statistics))Model with house, party and election type
mod_string2 = "
model {
## sampling
for (i in 1:N){
y[i] ~ dnorm(mu_party[PartyIndex[i]] + mu_house[HouseIndex[i]] + mu_elec[ElecIndex[i]] , invsigma2)
}
## priors
for (j in 1:J){
mu_party[j] ~ dnorm(mu, invtau2)
}
for (h in 1:H){
mu_house[h] ~ dnorm(mu, invtau2)
}
for (r in 1:R){
mu_elec[r] ~ dnorm(mu, invtau2)
}
invsigma2 ~ dgamma(1.5, 6)
sigma <- sqrt(pow(invsigma2, -1))
## hyperpriors
mu ~ dnorm(0, 4)
invtau2 ~ dgamma(3, 1)
tau <- sqrt(pow(invtau2, -1))
}"
set.seed(115)
params = c("mu_party", "mu_house","mu_elec", "mu", "invtau2", "invsigma2")
the_data2 <- list("y" = pt_bias_final$Bias,
"PartyIndex" = pt_bias_final$Partyiid,
"HouseIndex" = pt_bias_final$Companyiid,
"ElecIndex" = pt_bias_final$Eleciid,
"N" = length(pt_bias_final$Bias),
"J" = length(unique(pt_bias_final$Partyiid)),
"H" = length(unique(pt_bias_final$Companyiid)),
"R" = length(unique(pt_bias_final$Eleciid))
)mod2 = jags.model(textConnection(mod_string2), data=the_data2, n.chains=3)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 69
## Unobserved stochastic nodes: 18
## Total graph size: 354
##
## Initializing modelupdate(mod2, 1e3)
mod_sim2 = coda.samples(model=mod2,
variable.names=params,
n.iter=5e3)
mod_csim2 = as.mcmc(do.call(rbind, mod_sim2))
traceplot(mod_sim2, col = 1:3)
gelman.diag(mod_sim2)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## invsigma2 1 1
## invtau2 1 1
## mu 1 1
## mu_elec[1] 1 1
## mu_elec[2] 1 1
## mu_elec[3] 1 1
## mu_elec[4] 1 1
## mu_house[1] 1 1
## mu_house[2] 1 1
## mu_house[3] 1 1
## mu_house[4] 1 1
## mu_house[5] 1 1
## mu_party[1] 1 1
## mu_party[2] 1 1
## mu_party[3] 1 1
## mu_party[4] 1 1
## mu_party[5] 1 1
## mu_party[6] 1 1
##
## Multivariate psrf
##
## 1dic = dic.samples(mod2, n.iter=1e3)
(pm_params = colMeans(mod_csim2))## invsigma2 invtau2 mu mu_elec[1] mu_elec[2] mu_elec[3]
## 0.22210356 2.06050330 -0.30899280 -1.45596782 -0.25251530 0.67944498
## mu_elec[4] mu_house[1] mu_house[2] mu_house[3] mu_house[4] mu_house[5]
## -0.46195493 -0.22388919 -0.42778846 -0.25185377 0.06816059 -0.99148698
## mu_party[1] mu_party[2] mu_party[3] mu_party[4] mu_party[5] mu_party[6]
## 0.20019166 -0.35588102 -0.26834499 0.35541802 -1.25080284 -0.80848883summary(mod_sim2)##
## Iterations = 1001:6000
## Thinning interval = 1
## Number of chains = 3
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## invsigma2 0.22210 0.04162 0.0003399 0.0004504
## invtau2 2.06050 1.13126 0.0092367 0.0203188
## mu -0.30899 0.21647 0.0017675 0.0023309
## mu_elec[1] -1.45597 0.79879 0.0065221 0.0104836
## mu_elec[2] -0.25252 0.58265 0.0047573 0.0066304
## mu_elec[3] 0.67944 0.49682 0.0040565 0.0083672
## mu_elec[4] -0.46195 0.48444 0.0039554 0.0066320
## mu_house[1] -0.22389 0.54322 0.0044354 0.0060911
## mu_house[2] -0.42779 0.54539 0.0044531 0.0065991
## mu_house[3] -0.25185 0.56918 0.0046473 0.0061467
## mu_house[4] 0.06816 0.58492 0.0047759 0.0060925
## mu_house[5] -0.99149 0.46061 0.0037609 0.0065718
## mu_party[1] 0.20019 0.49779 0.0040644 0.0061182
## mu_party[2] -0.35588 0.68940 0.0056289 0.0064711
## mu_party[3] -0.26834 0.67890 0.0055432 0.0062635
## mu_party[4] 0.35542 0.70350 0.0057440 0.0081299
## mu_party[5] -1.25080 0.52259 0.0042669 0.0077507
## mu_party[6] -0.80849 0.47878 0.0039092 0.0063874
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## invsigma2 0.1486 0.1928 0.2191 0.24887 0.3109
## invtau2 0.6264 1.2565 1.8109 2.59124 4.9546
## mu -0.7307 -0.4505 -0.3108 -0.16725 0.1237
## mu_elec[1] -3.1928 -1.9394 -1.3959 -0.91002 -0.0566
## mu_elec[2] -1.3786 -0.6383 -0.2645 0.12766 0.9244
## mu_elec[3] -0.2374 0.3409 0.6546 0.99734 1.7289
## mu_elec[4] -1.4319 -0.7807 -0.4636 -0.15055 0.5043
## mu_house[1] -1.3102 -0.5769 -0.2185 0.12886 0.8571
## mu_house[2] -1.5252 -0.7821 -0.4180 -0.06726 0.6295
## mu_house[3] -1.3764 -0.6260 -0.2546 0.12432 0.8514
## mu_house[4] -1.0787 -0.3149 0.0578 0.45013 1.2319
## mu_house[5] -1.9251 -1.2882 -0.9765 -0.68052 -0.1219
## mu_party[1] -0.7609 -0.1285 0.1923 0.52150 1.2002
## mu_party[2] -1.7513 -0.7901 -0.3474 0.08937 0.9751
## mu_party[3] -1.6207 -0.7015 -0.2706 0.16700 1.0873
## mu_party[4] -0.9384 -0.1253 0.3136 0.80061 1.8541
## mu_party[5] -2.3336 -1.5840 -1.2307 -0.89727 -0.2795
## mu_party[6] -1.7774 -1.1211 -0.7978 -0.48500 0.1082bbb <-summary(mod_sim2)
# xtable(as.data.frame(bbb$statistics))Sensitivity analysis
For the sensitivity analysis we will assume very high values for our uncertain variance and nu_o
Model with just house and party
mod_string = "
model {
## sampling
for (i in 1:N){
y[i] ~ dnorm(mu_party[PartyIndex[i]] + mu_house[HouseIndex[i]] , invsigma2)
}
## priors
for (j in 1:J){
mu_party[j] ~ dnorm(mu, invtau2)
}
for (h in 1:H){
mu_house[h] ~ dnorm(mu, invtau2)
}
invsigma2 ~ dgamma(1.5, 15)
sigma <- sqrt(pow(invsigma2, -1))
## hyperpriors
mu ~ dnorm(0, 4)
invtau2 ~ dgamma(3, 1)
tau <- sqrt(pow(invtau2, -1))
}"
params = c("mu_party", "mu_house", "mu", "invtau2", "invsigma2")
the_data <- list("y" = pt_bias_final$Bias,
"PartyIndex" = pt_bias_final$Partyiid,
"HouseIndex" = pt_bias_final$Companyiid,
"N" = length(pt_bias_final$Bias),
"J" = length(unique(pt_bias_final$Partyiid)),
"H" = length(unique(pt_bias_final$Companyiid))
)
params = c("mu_party", "mu_house", "mu", "invtau2", "invsigma2")
mod = jags.model(textConnection(mod_string), data=the_data, n.chains=3)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 69
## Unobserved stochastic nodes: 14
## Total graph size: 260
##
## Initializing modelupdate(mod, 1e3)
mod_sim = coda.samples(model=mod,
variable.names=params,
n.iter=5e3)
mod_csim = as.mcmc(do.call(rbind, mod_sim))
(pm_params = colMeans(mod_csim))## invsigma2 invtau2 mu mu_house[1] mu_house[2] mu_house[3]
## 0.17612451 2.86514613 -0.29513064 -0.23349843 -0.36775627 -0.29418506
## mu_house[4] mu_house[5] mu_party[1] mu_party[2] mu_party[3] mu_party[4]
## 0.06359648 -0.90546652 0.15046624 -0.19755661 -0.13760190 0.03354543
## mu_party[5] mu_party[6]
## -1.05207361 -0.78294849mod_string = "
model {
## sampling
for (i in 1:N){
y[i] ~ dnorm(mu_party[PartyIndex[i]] + mu_house[HouseIndex[i]] , invsigma2)
}
## priors
for (j in 1:J){
mu_party[j] ~ dnorm(mu, invtau2)
}
for (h in 1:H){
mu_house[h] ~ dnorm(mu, invtau2)
}
invsigma2 ~ dgamma(1.5, 30)
sigma <- sqrt(pow(invsigma2, -1))
## hyperpriors
mu ~ dnorm(0, 4)
invtau2 ~ dgamma(3, 1)
tau <- sqrt(pow(invtau2, -1))
}"
params = c("mu_party", "mu_house", "mu", "invtau2", "invsigma2")
the_data <- list("y" = pt_bias_final$Bias,
"PartyIndex" = pt_bias_final$Partyiid,
"HouseIndex" = pt_bias_final$Companyiid,
"N" = length(pt_bias_final$Bias),
"J" = length(unique(pt_bias_final$Partyiid)),
"H" = length(unique(pt_bias_final$Companyiid))
)
params = c("mu_party", "mu_house", "mu", "invtau2", "invsigma2")
mod = jags.model(textConnection(mod_string), data=the_data, n.chains=3)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 69
## Unobserved stochastic nodes: 14
## Total graph size: 260
##
## Initializing modelupdate(mod, 1e3)
mod_sim = coda.samples(model=mod,
variable.names=params,
n.iter=5e3)
mod_csim = as.mcmc(do.call(rbind, mod_sim))
(pm_params = colMeans(mod_csim))## invsigma2 invtau2 mu mu_house[1] mu_house[2] mu_house[3]
## 0.16289180 2.86648265 -0.29582529 -0.23729588 -0.36226506 -0.30358782
## mu_house[4] mu_house[5] mu_party[1] mu_party[2] mu_party[3] mu_party[4]
## 0.04901296 -0.89755322 0.13634614 -0.20065987 -0.14169771 0.01642649
## mu_party[5] mu_party[6]
## -1.03181489 -0.76678259Model with house, party and election type
mod_string2 = "
model {
## sampling
for (i in 1:N){
y[i] ~ dnorm(mu_party[PartyIndex[i]] + mu_house[HouseIndex[i]] + mu_elec[ElecIndex[i]] , invsigma2)
}
## priors
for (j in 1:J){
mu_party[j] ~ dnorm(mu, invtau2)
}
for (h in 1:H){
mu_house[h] ~ dnorm(mu, invtau2)
}
for (r in 1:R){
mu_elec[r] ~ dnorm(mu, invtau2)
}
invsigma2 ~ dgamma(1.5, 15)
sigma <- sqrt(pow(invsigma2, -1))
## hyperpriors
mu ~ dnorm(0, 4)
invtau2 ~ dgamma(3, 1)
tau <- sqrt(pow(invtau2, -1))
}"
set.seed(115)
params = c("mu_party", "mu_house","mu_elec", "mu", "invtau2", "invsigma2")
the_data2 <- list("y" = pt_bias_final$Bias,
"PartyIndex" = pt_bias_final$Partyiid,
"HouseIndex" = pt_bias_final$Companyiid,
"ElecIndex" = pt_bias_final$Eleciid,
"N" = length(pt_bias_final$Bias),
"J" = length(unique(pt_bias_final$Partyiid)),
"H" = length(unique(pt_bias_final$Companyiid)),
"R" = length(unique(pt_bias_final$Eleciid))
)
mod2 = jags.model(textConnection(mod_string2), data=the_data2, n.chains=3)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 69
## Unobserved stochastic nodes: 18
## Total graph size: 354
##
## Initializing modelupdate(mod2, 1e3)
mod_sim2 = coda.samples(model=mod2,
variable.names=params,
n.iter=5e3)
mod_csim2 = as.mcmc(do.call(rbind, mod_sim2))
traceplot(mod_sim2, col = 1:3)
gelman.diag(mod_sim2)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## invsigma2 1 1
## invtau2 1 1
## mu 1 1
## mu_elec[1] 1 1
## mu_elec[2] 1 1
## mu_elec[3] 1 1
## mu_elec[4] 1 1
## mu_house[1] 1 1
## mu_house[2] 1 1
## mu_house[3] 1 1
## mu_house[4] 1 1
## mu_house[5] 1 1
## mu_party[1] 1 1
## mu_party[2] 1 1
## mu_party[3] 1 1
## mu_party[4] 1 1
## mu_party[5] 1 1
## mu_party[6] 1 1
##
## Multivariate psrf
##
## 1dic = dic.samples(mod2, n.iter=1e3)
(pm_params = colMeans(mod_csim2))## invsigma2 invtau2 mu mu_elec[1] mu_elec[2] mu_elec[3]
## 0.20759134 2.16131742 -0.31059633 -1.38158267 -0.26186068 0.63429654
## mu_elec[4] mu_house[1] mu_house[2] mu_house[3] mu_house[4] mu_house[5]
## -0.47479163 -0.21559704 -0.41189908 -0.24463000 0.04654078 -0.96681012
## mu_party[1] mu_party[2] mu_party[3] mu_party[4] mu_party[5] mu_party[6]
## 0.18443221 -0.33279075 -0.26750023 0.30600927 -1.21179744 -0.78663631mod_string2 = "
model {
## sampling
for (i in 1:N){
y[i] ~ dnorm(mu_party[PartyIndex[i]] + mu_house[HouseIndex[i]] + mu_elec[ElecIndex[i]] , invsigma2)
}
## priors
for (j in 1:J){
mu_party[j] ~ dnorm(mu, invtau2)
}
for (h in 1:H){
mu_house[h] ~ dnorm(mu, invtau2)
}
for (r in 1:R){
mu_elec[r] ~ dnorm(mu, invtau2)
}
invsigma2 ~ dgamma(1.5, 30)
sigma <- sqrt(pow(invsigma2, -1))
## hyperpriors
mu ~ dnorm(0, 4)
invtau2 ~ dgamma(3, 1)
tau <- sqrt(pow(invtau2, -1))
}"
set.seed(115)
params = c("mu_party", "mu_house","mu_elec", "mu", "invtau2", "invsigma2")
the_data2 <- list("y" = pt_bias_final$Bias,
"PartyIndex" = pt_bias_final$Partyiid,
"HouseIndex" = pt_bias_final$Companyiid,
"ElecIndex" = pt_bias_final$Eleciid,
"N" = length(pt_bias_final$Bias),
"J" = length(unique(pt_bias_final$Partyiid)),
"H" = length(unique(pt_bias_final$Companyiid)),
"R" = length(unique(pt_bias_final$Eleciid))
)
mod2 = jags.model(textConnection(mod_string2), data=the_data2, n.chains=3)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 69
## Unobserved stochastic nodes: 18
## Total graph size: 354
##
## Initializing modelupdate(mod2, 1e3)
mod_sim2 = coda.samples(model=mod2,
variable.names=params,
n.iter=5e3)
mod_csim2 = as.mcmc(do.call(rbind, mod_sim2))
traceplot(mod_sim2, col = 1:3)
gelman.diag(mod_sim2)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## invsigma2 1 1
## invtau2 1 1
## mu 1 1
## mu_elec[1] 1 1
## mu_elec[2] 1 1
## mu_elec[3] 1 1
## mu_elec[4] 1 1
## mu_house[1] 1 1
## mu_house[2] 1 1
## mu_house[3] 1 1
## mu_house[4] 1 1
## mu_house[5] 1 1
## mu_party[1] 1 1
## mu_party[2] 1 1
## mu_party[3] 1 1
## mu_party[4] 1 1
## mu_party[5] 1 1
## mu_party[6] 1 1
##
## Multivariate psrf
##
## 1dic = dic.samples(mod2, n.iter=1e3)
(pm_params = colMeans(mod_csim2))## invsigma2 invtau2 mu mu_elec[1] mu_elec[2] mu_elec[3]
## 0.18720907 2.30782506 -0.30919797 -1.26847525 -0.27761760 0.57370895
## mu_elec[4] mu_house[1] mu_house[2] mu_house[3] mu_house[4] mu_house[5]
## -0.48394649 -0.20971232 -0.39773000 -0.24919981 0.03576236 -0.93825866
## mu_party[1] mu_party[2] mu_party[3] mu_party[4] mu_party[5] mu_party[6]
## 0.16292985 -0.32405373 -0.25471263 0.25056884 -1.15488402 -0.76468216